function model = ProjLoadingsExtraControls(model)

% Initial call to ProjectExtraControls to get the dimension of yExtra
[yExtra0_cu,namesExtra,yExtraSS] = ProjectExtraControls(zeros(model.ny,1),zeros(model.nx,1),model);
nyExtra    = length(namesExtra);
xGrid      = model.xGrid;
nGrid      = size(xGrid,1);
yExtra     = zeros(nGrid,nyExtra);
for i=1:nGrid
    x_cu = model.xGrid(i,:)';
    if model.appOrder == 1
        y_cu   = model.g0 + model.gx*x_cu;    
    elseif model.appOrder == 2
        AA_cu  = kron3(x_cu,x_cu);
        y_cu   = model.g0 + model.gx*x_cu + model.Gxx*AA_cu;    
    elseif model.appOrder == 3
        AA_cu  = kron3(x_cu,x_cu);
        AAA_cu = kron3(AA_cu,x_cu);
        y_cu   = model.g0 + model.gx*x_cu + model.Gxx*AA_cu + model.Gxxx*AAA_cu;    
    end
    tmpy   = ProjectExtraControls(y_cu,x_cu,model);
    yExtra(i,:) = struct2array(tmpy);
end
nx = model.nx;
% The regression - done without scaling of the states (not needed)
[theta,resUnitFree] = OLSforProj(yExtra,model.xGrid,model.appOrder,model.scaleX_in_g*0+1);
%% Unfolding theta
g0 = theta(:,1);

% First order terms
gx  = theta(:,1+1:1+nx);

% Second order terms
gxx    = zeros(nyExtra,nx,nx);
if  model.appOrder >= 2
    idx    = 1+nx;
    for i=1:nx
        for j=i:nx
            idx = idx + 1;
            gxx(:,i,j) = theta(:,idx);
            gxx(:,j,i) = gxx(:,i,j);
        end
    end
end
% Check
%gFit = zeros(nyExtra,nGrid);
%Gxx = reshape(gxx,nyExtra,nx^2);
%for i=1:nGrid
%    gFit(:,i) = g0 + gx*model.xGrid(i,:)' + Gxx*kron(model.xGrid(i,:)',model.xGrid(i,:)');
%end
%check2nd = [gFit-Yfit]

% Third order terms
gxxx  = zeros(nyExtra,nx,nx,nx);
if  model.appOrder >= 3
    for alfa1=1:nx
        for alfa2=alfa1:nx
            for alfa3=alfa2:nx
                idx = idx + 1;
                gxxx(:,alfa1,alfa2,alfa3) = theta(:,idx);
                
                % Using symmetry for alfa1 and alfa2
                if alfa1 == alfa2 && alfa2 ~= alfa3 %alfa1==alfa2~=alfa3
                    %gxxx(:,alfa1,alfa1,alfa3)= gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa1,alfa3,alfa1) = gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa3,alfa1,alfa1) = gxxx(:,alfa1,alfa2,alfa3);
                end
                % Using symmetry for alfa2 and alfa3
                if alfa1 ~= alfa2 && alfa2 == alfa3  %alfa1~=alfa2==alfa3
                    %gxxx(:,alfa1,alfa2,alfa2)= gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa2,alfa1,alfa2) = gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa2,alfa2,alfa1) = gxxx(:,alfa1,alfa2,alfa3);
                end
                % Using symmetry for alfa1,alfa2, and alfa3
                if alfa1 ~= alfa2 && alfa1 ~= alfa3 &&  alfa2 ~= alfa3 %alfa1~=alfa2~=alfa3
                    %gxxx(:,alfa1,alfa2,alfa3) = gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa1,alfa3,alfa2) = gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa3,alfa1,alfa2) = gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa3,alfa2,alfa1) = gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa2,alfa3,alfa1) = gxxx(:,alfa1,alfa2,alfa3);
                    gxxx(:,alfa2,alfa1,alfa3) = gxxx(:,alfa1,alfa2,alfa3);
                end
            end
        end
    end
end
% Check
%gFit = zeros(nyExtra,nGrid);
%Gxx = reshape(gxx,nyExtra,nx^2);
%Gxxx = reshape(gxxx,nyExtra,nx^3);
%for i=1:nGrid
%    AA = kron(model.xGrid(i,:)',model.xGrid(i,:)');
%    gFit(:,i) = g0 + gx*model.xGrid(i,:)' + Gxx*AA + Gxxx*kron(model.xGrid(i,:)',AA);
%end
%check3rd = max(max(abs([gFit-Yfit])))
            
%% Saving the results
model.labely = [model.labely namesExtra];
model.ny     = length(model.labely);
model.gSS    = [model.gSS;struct2array(yExtraSS)'];
model.g0     = [model.g0;g0];
model.gx     = [model.gx;gx];
model.gxx    = [model.gxx;gxx];
model.gxxx   = [model.gxxx;gxxx];
model.Gxx    = reshape(model.gxx,model.ny,nx^2);
model.Gxxx   = reshape(model.gxxx,model.ny,nx^3);
model.residualsExtra = resUnitFree;

end
